gc()
rm(list = ls())
library(yyeasy)
library(extrafont)
yyload(ggalluvial, RColorBrewer, pheatmap, cowplot)

group2 <- yyread("0_data/all_data.xlsx", excel = T, sheet = "编号2")
rownames(group2) <- group2$id3

# gene_tpm <- yyread("0_data_meta/gene.tpm.tsv", rownames = T)
# gene_count <- yyread("0_data_meta/gene.count.tsv", rownames = T)
sp_count <- yyread("0_data_meta/sp.count.tsv", rownames = T)
# yywrite(
#   select(sp_count , -taxonomy), "0_data_meta/sp.count_no_tax.tsv", row.names = T
# )
sp_asv <- sp_count[, 1:9]
sp_tax <- sp_count[, 10, drop = F] %>%
  separate(1, c("k", "p", "c", "o", "f", "g", "s"), sep = "; ")
sp_phy <- phyloseq(
  otu_table(sp_asv, T),
  tax_table(as.matrix(sp_tax))
)
### 折叠----------------------------------
identical(rownames(otu_table(sp_phy)), rownames(tax_table(sp_phy)))
sp_phy_bind <- bind_cols(otu_table(sp_phy), tax_table(sp_phy))
sp_ge <- sp_phy_bind %>%
  select(-s) %>%
  group_by(k, p, c, o, f, g) %>%
  summarise_all(sum) %>%
  ungroup() %>%
  select(-1:-4)
rep_id <- sp_ge$g == "g__"
sp_ge$g[rep_id] <- paste0(sp_ge$f[rep_id], "; ", sp_ge$g[rep_id])
sp_ge <- select(sp_ge, -f) %>% column_to_rownames("g")




### p-----------------
sp_ph <- sp_phy_bind %>%
  select(-c(c, o, f, g, s)) %>%
  group_by(k, p) %>%
  summarise_all(sum) %>%
  ungroup() %>%
  select(-1) %>%
  column_to_rownames("p")
## k ----------------------
### p-----------------
sp_ki <- sp_phy_bind %>%
  select(-c(p, c, o, f, g, s)) %>%
  group_by(k) %>%
  summarise_all(sum) %>%
  column_to_rownames("k")

rownames(sp_ph) <- str_sub(rownames(sp_ph), 4, -1)
rownames(sp_ge) <- str_sub(rownames(sp_ge), 4, -1)
rownames(sp_ki) <- str_sub(rownames(sp_ki), 4, -1)

dat <- sp_ph
N <- 15
##
top_name <- rowSums(dat) %>%
  sort(T) %>%
  names()
if (length(top_name) > N) {
  top1 <- dat[top_name[1:N], ]
  top2 <- rbind(
    top1,
    Others = colSums(dat) - colSums(top1)
  )
  top_total <- decostand(top2, "total", 2) %>% rownames_to_column("Taxon")
} else {
  top_total <- dat[top_name, ] %>%
    decostand("total", 2) %>%
    rownames_to_column("Taxon")
}
# tt_asv <- (column_to_rownames(top_total, "Taxon")*1000) %>% select(contains("玉米"), contains("水稻"), contains("空白"))
# tt_ge  <- (column_to_rownames(top_total, "Taxon")*1000) %>% select(contains("玉米"), contains("水稻"), contains("空白"))
# # tt_py  <- (column_to_rownames(top_total, "Taxon")*1000) %>% select(contains("玉米"), contains("水稻"), contains("空白"))
# yywrite(tt_asv, "tt_sp.tsv", row.names = T)
# yywrite(tt_ge, "tt_ge.tsv", row.names = T)
# yywrite(tt_py, "tt_py.tsv", row.names = T)

## 生成绘图数据
datBar <- pivot_longer(top_total, -Taxon, names_to = "sample", values_to = "count")
datBar$Taxon <- factor(datBar$Taxon, levels = rev(top_total$Taxon))
datBar$sample <- factor(datBar$sample, levels = group2$id3)
datBar2 <- filter(datBar, Taxon != "Others")

Colors <- c()
for (color in c("Set1", "Set2", "Set3", "Paired", "Dark2")) {
  Colors <- c(Colors, brewer.pal(8, color))
}
tmp_len <- length(unique(datBar$Taxon))

## 设置 Taxon 的因子顺序并绘制barplot
### phylum -----------------------
(p1 <- ggplot(datBar, aes(x = sample, y = count, alluvium = Taxon, stratum = Taxon)) +
  geom_col(aes(fill = Taxon), width = 0.6) +
  geom_alluvium(aes(fill = Taxon), alpha = .5) +
  labs(x = NULL, y = "Sequence Percent", fill = "Species") +
  scale_fill_manual(values = Colors[][1:tmp_len]) +
  scale_y_continuous(expand = c(0, 0)) +
  # coord_cartesian(ylim = c(0.9, 1)) +
  theme_classic2(10.5) +
  theme(
    axis.title = element_text(family = "serif"),
    axis.text.y = element_text(family = "serif"),
    axis.text.x = element_text(family = "SimSun", angle = 45, vjust = 0.5),
    legend.text = element_text(family = "serif", size = 10, face = "italic"),
    legend.key.size = unit(0.18, "inches"),
    legend.position = "none"
  ))
(p2 <- ggplot(datBar, aes(x = sample, y = count, alluvium = Taxon, stratum = Taxon)) +
  geom_col(aes(fill = Taxon), width = 0.6) +
  geom_alluvium(aes(fill = Taxon), alpha = .5) +
  labs(x = NULL, y = "Sequence Percent", fill = "Phylum") +
  scale_fill_manual(values = Colors[][1:tmp_len]) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(ylim = c(0.9, 1)) +
  theme_classic2(10.5) +
  theme(
    axis.title = element_text(family = "serif"),
    axis.title.y = element_blank(),
    axis.text.y = element_text(family = "serif"),
    axis.text.x = element_text(family = "SimSun", angle = 45, vjust = 0.5),
    legend.text = element_text(family = "serif", size = 10),
    legend.key.size = unit(0.14, "inches")
  ))
p3 <-
  cowplot::plot_grid(
    p1,
    p2 + theme(legend.position = "none"),
    cowplot::get_legend(p2),
    rel_widths = c(1, 1, .48),
    nrow = 1,
    labels = c("a", "b", ""),
    label_fontfamily = "serif",
    label_fontface = "plain",
    label_y = .98,
    label_x = -0.01
  )
p3


yysave(p3, "sp_ph2.png", 18, 10, dpi = 1200)
yysave(p3, "sp_ph2.svg", 18, 10)

## --------genus------------------------------
(p4 <- ggplot(datBar, aes(x = sample, y = count, alluvium = Taxon, stratum = Taxon)) +
  geom_col(aes(fill = Taxon), width = 0.6) +
  geom_alluvium(aes(fill = Taxon), alpha = .5) +
  labs(x = NULL, y = "Sequence Percent", fill = "Genus") +
  scale_fill_manual(values = Colors[][1:tmp_len]) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(ylim = c(0, 0.7)) +
  theme_classic2(10.5) +
  theme(
    axis.title = element_text(family = "serif"),
    axis.text.y = element_text(family = "serif"),
    axis.text.x = element_text(family = "SimSun", angle = 45, vjust = 0.5),
    legend.text = element_text(family = "serif", size = 10, face = "italic"),
    legend.key.size = unit(0.18, "inches")
  ))
yysave(p4, "sp_ge2.png", 15, 10, dpi = 1200)
yysave(p4, "sp_ge2.svg", 15, 10)



## kingdom
datBar$Taxon <- factor(datBar$Taxon, levels = (top_total$Taxon))
(p5 <- ggplot(datBar, aes(x = sample, y = count, alluvium = Taxon, stratum = Taxon)) +
  geom_col(aes(fill = Taxon), width = 0.6) +
  geom_alluvium(aes(fill = Taxon), alpha = .5) +
  labs(x = NULL, y = "Sequence Percent", fill = "Kingdom") +
  scale_fill_manual(values = Colors[][1:tmp_len]) +
  scale_y_continuous(expand = c(0, 0)) +
  coord_cartesian(ylim = c(0, 0.04)) +
  theme_classic2(10.5) +
  theme(
    axis.title = element_text(family = "serif"),
    axis.text.y = element_text(family = "serif"),
    axis.text.x = element_text(family = "SimSun", angle = 45, vjust = 0.5),
    legend.text = element_text(family = "serif", size = 10),
    legend.key.size = unit(0.18, "inches")
  ))

yysave(p5, "sp_ki4.png", 15, 10, dpi = 1200)
yysave(p5, "sp_ki4.svg", 15, 10)

###          alpha beta dist-------------------------------------
# fy_df <- rarefy_otu(sp_phy)
# save(fy_df, file = "0_data/fy_df.Rdata")
x_min <- group_by(fy_df, id) %>%
  summarise(across(Depth, max)) %>%
  .[, "Depth"] %>%
  min()
p_rary <- ggplot(fy_df, aes(Depth, Observed, group = id, color = id)) +
  geom_vline(xintercept = x_min, linetype = "dashed", color = "grey60") +
  stat_summary(fun = "mean", geom = "smooth", size = .7) +
  stat_summary(
    fun.data = "mean_cl_normal", geom = "errorbar",
    size = .5, aes(width = max(Depth) / 30)
  ) +
  labs(y = "Observed species", x = "Sequencing depth") +
  theme_bw2(10.5) +
  theme(
    legend.justification = c(1, 0),
    legend.position = c(1, 0),
    axis.title = element_text(family = "serif"),
    axis.text.y = element_text(family = "serif"),
    axis.text.x = element_text(family = "serif", angle = 0, vjust = 1),
    legend.text = element_text(family = "serif", size = 10),
    legend.key.size = unit(0.18, "inches")
  )
yysave(p_rary, "p_rare.svg", 12, 10)
yysave(p_rary, "p_rare.png", 12, 10, dpi = 1200)


sp_alpha <- estimate_richness(sp_phy)
sp_alpha2 <- alpha_index(sp_asv) %>% select(-4, -5)
yywrite(sp_alpha2, "tmp.tsv", row.names = T)
sp_dist2 <- beta_dist(sp_asv)

p_beta <- pcoa_plot(
  sp_dist2$Bray_Curtis,
  data.frame(
    colnames(sp_asv),
    str_sub(colnames(sp_asv), 1, 4)
  ),
  font = "SimSun",
  font_size = 10.5
)$p +
  aes(color = group) +
  theme_bw2(10.5) +
  theme(
    axis.title = element_text(family = "serif"),
    axis.text = element_text(family = "serif"),
    legend.text = element_text(family = "serif", size = 10),
    legend.position = "none"
  )

anosim(sp_dist2$Bray_Curtis, str_sub(colnames(sp_asv), 1, 4))
adonis(sp_dist2$Bray_Curtis ~ str_sub(colnames(sp_asv), 1, 4))
## anosim: R = 1, P = 0.008
## adonis: R2 = 0.778, P =  0.003


yysave(p_beta, "sp_beta.png", 12, 10, dpi = 1200)
yysave(p_beta, "sp_beta.svg", 12, 10)

### pheat
tmp_cor <- psych::corr.test(sp_asv)
tmp_dist <- as.dist(1 - tmp_cor$r)


pp <- pheat(tmp_cor, , "holm", fontfamily = "SimSun", fontsize = 10, treeheight_row = 30, treeheight_col = 30, angle_col = 45, clustering_distance_rows = tmp_dist, clustering_distance_cols = tmp_dist)

yysave(pp, "pheat_sp.svg", 12, 11)
yysave(pp, "pheat_sp.png", 12, 11, dpi = 1200)

### 薇恩图